Constraints on local primordial non-Gaussianity from large scale structure 
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Recent work has shown that the local non-Gaussianity parameter /nl induces a scale-dependent 
bias, whose amplitude is growing with scale. Here we first rederive this result within the context of 
peak-background split formalism and show that it only depends on the assumption of universality 
of mass function, assuming halo bias only depends on mass. We then use extended Press-Schechter 
formalism to argue that this assumption may be violated and the scale dependent bias will depend 
on other properties, such as merging history of halos. In particular, in the limit of recent mergers 
we find the effect is suppressed. Next we use these predictions in conjunction with a compendium 
of large scale data to put a limit on the value of /nl- When combining all data assuming that halo 
occupation depends only on halo mass, we get a limit of —29 (—65) < /nl < -|-70 (-1-93) at 95% 
(99.7%) confidence. While we use a wide range of datasets, our combined result is dominated by the 
signal from the SDSS photometric quasar sample. If the latter are modelled as recent mergers then 
the limits weaken to —31 (—96) < /nl < +70 (+96). These limits are comparable to the strongest 
current limits from the WMAP 5 year analysis, with no evidence of a positive signal in /nl. While 
the method needs to be thoroughly tested against large scale structure simulations with realistic 
quasar and galaxy formation models, our results indicate that this is a competitive method relative 
to CMB and should be further pursued both observationally and theoretically. 

PACS numbers; 98.80.Jk, 98.80.Cq 



I. INTRODUCTION 

The origin of structure formation in the universe is 
one of the most hotly debated topics in current cosmol- 
ogy research. The standard paradigm is that of inflation 
[J S B which has been tremendously successful in 
describinga very large number of very distinct data-sets 
(see e.g. Q). Inflationary models generically predict a 
flat universe and nearly scale-invariant spectrum of ini- 
tial fluctuations [1, 0, S, H, US] , both of which seem to be 
confirmed by observations. Consequently, a lot of effort is 
being put into constraining observables that might actu- 
ally distinguish between different models of inflation. At 
the moment, this is done as a multi-pronged effort: first, 
measurement of the primordial power spectrum gives a 
direct measure of the inflationary potential shape and 
inflationary models differ on the actual slope, some pre- 
dicting red and some blue spectrum. Moreover, inflation 
predicts that the primordial slope should only be chang- 
ing with scale very slowly and any deviation from this 
prediction would be a surprise in need of an explana- 
tion. Second, a detection of _B-mode polarization in the 
cosmic microwave background, if interpreted as gravita- 
tional waves from the early Universe, will effectively de- 
termine the energy scale of inflation and rule out a major 
class of inflationary models that predict inflation occurs 
at a low energy scale 111 
as ekpyrotic models (ill 



Alternatives to inflation, such 
1^ , [T3 | , differ from inflation- 



ary predictions in that the expected gravitational wave 
signal in CMB is always negligible. Thus, they can be 
falsifled if primordial gravitational waves are detected. 
Third, multifield models could generate isocurvature per- 
turbations; while now ruled out as the main mode of 
structure formation, these could be present at a subdom- 
inant level and if detected would rule out the simplest 
models of inflation [H, [11, [13, [13 . 

A fourth direction, and the one we focus on in this 
paper, is non-Gaussianity in initial conditions. Stan- 
dard single field inflation predicts that the departures 
from Gaussianity are very small and not accessible to 
the current observational constraints. Most of the models 
predict that non-Gaussianity is of the local type, mean- 
ing that it depends on the local value of the potential 
only. A standard parameterization of the primordial 
non-Gaussianity is the so-called scale-independent /nl 
parameterization, in which one includes a quadratic cor- 
rection to the potential [l^, H^l : 



(1) 



where cf) is the primordial potential assumed to be a Gaus- 
sian random field and /nl describes the amplitude of the 
correction. A typical value of /nl for standard slow roll 
inflation is of the order of slow roll parameter and thus 
of order 
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2l| . but this is likely to be swamped by 



the contribution from nonlinear transformation between 
the primordial field fluctuation (assumed to be Gaussian 
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if it started from the pure Bunch-Davies vacuum) and 
the observable (such as CMB temperature fluctuation), 
which generically gives /nl of order unity (see e.g. [2^1). 
Models where is significantly higher include multi- 
field inflation [2^ [2^ as well as models where non- 
Gaussianity arises during reheating [iH, H^l or preheating 
[28l . [29l . [30l . [sH . Non-slow roll inflation models may also 
lead to a significant non-Gaussianity, but are constrained 
because they may not lead to enough inflation in the first 
place. Note that since (j> ~ iO~^ even /nl ~ 100, compa- 
rable to the present limits, generates non-Gaussian sig- 
natures only at a 10"'^ level, so the non-Gaussian signal 
one is searching for is very small. Overall, any detection 
of /nl above unity would be a major surprise in need of 
an explanation within the inflationary paradigm. 

Very recently, non-Gaussianity in ekpyrotic models has 
also been studied with the results suggesting that non- 
Gaussianity in these models is generically lar ge 13211 and 
often correlated with the spectral slope |33l . |3J|, in 
the sense that the redder the spectrum the higher non- 
Gaussianity one may expect. Thus, non-Gaussianity is 
emerging as one of the strongest discriminators among 
the models attempting to explain the origins of structure 
in the universe. 

Traditionally, the cleanest method for detecting the 
non-Gaussianity has been to measure the bispectrum or 
3-point function of the Cosmic Microwave Background 
(CMB). The initial 3-year Wilkinson Anisotropy Probe 
(WMAP) result gave a hmit on /nl of -54 < /nl < 134 
[351] (all limits reported at 95% confidence limit) from 
the bispectrum of the WMAP data at ^ < 400. This has 
been improved subsequently to —36 < /nl < 100 by the 
3 year analysis [36]. However, Yadav & Wandelt recently 
claimed a detection of /nl > at 99.5% significance, with 
2(7 range 27 < /nl < 147 [? ]. This result is surprising 
as, taken at its face value, it implies a very non-standard 
inflation or something else all together. Interestingly, a 
similar result has been obtained by Jeong & Smoot using 
the one-point distribution function of the CMB [sl] . The 
WMAP 5-year results also prefer positive /nl from their 
bispectrum analysis, although zero /nl is within its 2- 
sigma significance [a], namely — 9 < /nl < 111. 

An alternative method that has been applied to CMB 
are Minkowski functionals [35i]. WMAP 3-year analysis 
via this technique puts the limit at — 70 < /nl < 91 [39j . 
while the recent WMAP 5-year analysis gives —178 < 
/nl < 64 from Minkowski functionals, which is about a 
factor of two larger error than from bispectrum analysis 
i. 

In the near future, the Planck satellite should improve 
these numbers significantly and can in principle push 
to tT(/NL) ^ 7 [40|, while more speculatively, the pre- 
reionization H i 21 cm transition may offer an unprece- 
dented access to the 3-dimensional distribution of linear 
modes at high redshift and may bring us into regime 
('■(/nl) ^ 1 where a detection is expected 141). 

A different direction to probe non-Gaussianity is to try 
and determine what observational signatures it leaves in 



the large scale structure (LSS) of the universe. The main 
problem is that non-linearities add their own phase cor- 
relations between Fourier modes that can very quickly 
swamp the primordial signal. Historically, the focus was 
on the mass function of very massive virialised structures 
m, m, m, m Si] . The motivation was the notion that 
very massive virialised objects correspond to very rare 
peaks in the initial density field and therefore their num- 
ber density should be an exponentially sensitive probe 
of those peaks at the high mass end, allowing a unique 
probe of the primordial peak structure. While results 
generally agree with this picture, the observational task 
is made very difficult by the low number statistics of such 
objects, uncertainties in the mass-observable relation and 
its scatter, and selection effects. 

A different method has been recently proposed by 
Dalai et al. [13]. By extending the classical calculation 
for calculating the clustering of rare peaks in a Gaus- 
sian field [11] to the /NL-type non-Gaussianity, they have 
shown that clustering of rare peaks exhibits a very dis- 
tinct scale-dependent bias on the largest scales. The an- 
alytical result has been tested using iV-body simulations, 
which confirm this basic picture. 

The purpose of this paper is two-fold, first to provide 
a better theoretical understanding of the effect and the 
range of its applicability, and second, to apply it to the 
real data. We begin in Section [H] by providing a new, 
more general, derivation of the nonlinear bias induced by 
non-Gaussianity, highlighting more clearly its underlying 
assumptions. We then extend the basic derivation using 
the extended Press-Schechter formalism and show that 
for certain classes of halos, such as those that have un- 
dergone a recent merger, the results may be substantially 
modified. 

We then use this formalism and apply it to a wide se- 
lection of publicly available large scale structure data. In 
Sections IIIIHIVI we discuss the data, methodology, main 
results and systematic issues, including application of 
Section |TT] to the derived observational constraints. In 
Section |V] we discuss the results and present some direc- 
tions for the future. 



II. THEORY 

In this section we provide theoretical derivations of the 
large scale bias induced by non-Gaussianity of the local 
type. We first derive an expression that depends only 
on the halo mass function and halo bias, using the fact 
that /nl causes a local re-scaling of the amplitude erg. 
This derivation is more general than that of Dalai et al. 
[13], since it is not tied to the spherical collapse model. 
In particular, we show that any universal mass function, 
such as the Sheth & Tormen mass function or Press- 
Schechter mass functio n ISOl . leads to the equation first 
derived in Dalai et al. [43]. We then extend the deriva- 
tion to the extended Press-Schechter (ePS) type of anal- 
ysis and derive the effect of halo merger bias on Ab. Pi- 
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nally we comment on the accuracy of the ePS prediction 
and compare it to previously pubUshed A^-body results. 



Local non-Gaussianity in peak-background 
formalism 



Large-scale bias of haloes is usually treated in the con- 
text of the peak-background split [5l| . One can split the 
density field into a long-wavelength piece Si and a short- 
wavelength piece Ss as in 



(2) 



The local Lagrangian number density of haloes n(x) 
(i.e. number density of haloes per unit halo mass) at 
position X can then be written as a function of the lo- 
cal value of the long- wavelength perturbation Si{x.) and 
the statistics of the short- wavelength fluctuations Ps{ks). 
The sufficiently averaged local density of halos follows the 
large scale matter perturbations 



n(x) =n{l + BlSi) 
and so the Lagrangian bias is then 



(3) 



(4) 



For Eulerian space bias one needs to add the Eule- 
rian space clustering, so the total or Eulerian bias is 
6 = 6l + 1. This argument leads to a generically scale- 
independent bias at sufficiently large scales. The specific 
function 6(M) is obtained by constructing a specific func- 
tion n[6i {yi), Ps{ks); M], generally fit to simulations, and 
then differentiating it. 

The non-Gaussian case is complicated by the fact that 
large and small-scale density fluctuations are no longer 
independent. Instead, in the /nl prescription, one may 
separate long- and short-wavelength Gaussian potential 
fluctuations, 



(5) 



which are independent. Inserting into Equation ([T]) we 
can then re-map these into the non-Gaussian potential 
fluctuations, 

* = fNLcjjf + (1 + 2/nl</'/)?!'s + InlcI^I + const. (6) 

We can then convert this to a density field using the 
expression 6i{k) = a{k)^{k), with 



a(k) 



2c^PT{k)D{z) 



3a 



(7) 



Here T(fc) is the transfer function, c speed of light, D{z) 
the linear growth factor normalised to be (1 + z)^^ in the 
matter domination, the matter density today and Hq 
the Hubble parameter today. The operator a{k) makes 



it non-local on scales of 100 Mpc, so this can also be 
thought of as a convolution operator in real space. 

For long-wavelength modes of the density field, one 
may write 

6i{k) = a{k)Uky, (8) 

the remaining terms in Equation ([6]) are either much 
smaller (/nl^^ ); have only short- wavelength pieces [(1-1- 
2/nl(/);)0s], or simply add a small white noise contribu- 
tion on large scales (/nl'/'s)- 

Within a region of given large-scale over-density di and 
potential (pi, the short- wavelength modes of the density 
field are: 

S^^a [(1 + 2hM<Ps + hL(t>l] ■ 

This is a special case of 

d,=a [Xi^s + ^202] , 

where Xi = 1 + 2f^i^4>i and X2 = /nl- 

In the non-Gaussian case, the local number density of 
haloes of mass M is a function of not just 61, but also Xi 
and X2: n[6i, Xi, X2; Ps{ks)- M]. The halo bias is then 



(9) 



(10) 



bL{M,k) = 



dr 



2/1 



_d6i{^) • ^■''''''dSiik) dXi 

where the derivative is taken at the mean value Xi = 1. 
(There is no X2 term since X2 is not spatially variable.) 
The first term here is the usual Gaussian bias, which has 
no dependence on k. 

Equation IjlOp shows that the effect on non-Gaussianity 
is a local rescaling of amplitude of (small scale) matter 
fluctuations. To keep the cosmologist's intuition we write 
this in terms of cts: 



d(j)i{k) dn 



(11) 



4°^-'(x) =(78Xi(x), 



(12) 



so S(tI°'"'^ 
PT|) as 



asSXi. This allows us to rewrite Equation 



d4>i{k) dlnn 
dSi{k) d\nai°'^^' 



(13) 



In principle there is an additional change in the bias 
because the mean density h contains terms of order /nl, 
which arise from (i) the dependence of n on X2 and (ii) 
the cross-correlation of Si and Xi. This correction is 
scale-independent and so cause no problem if one is flt- 
ting the bias to large-scale structure data, as we do here. 

Substituting in d(j>i{k)/ddi{k) = a^^{k) and dropping 
the local label, we flnd: 



Ab{M, k) 



3a 



dlnn 



c?BT{k)D{zy^^'d\na8 



(14) 



This formula is extremely useful because it applies to 
the bias of any type of object and is expressible entirely 
in terms of quantities in Gaussian cosmologies, which 
have received enormous attention from A'^-body simula- 
tors. Within the peak-background split model, the task 
of performing non-Gaussian calculations is thus reduced 
to an ensemble of Gaussian simulations with varying am- 
plitude of matter fluctuations. 
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B. Application to universal mass functions 



C. Halo merger bias 



We now apply Equation to halo abundance models 
with a universal mass function. Universal mass functions 
are those that depend only significance i^(M), i.e. 



n{M) = n{M, v) = M-^vf{v) 



d In M ' 



(15) 



where we define v = 5l/(P{M) and f{y) is the fraction 
of mass that collapses into haloes of significance between 
V and V + dv. Here 5c — 1.686 denotes the spherical 
collapse linear over-density and a{M) is the variance of 
the density field smoothed with a top-hat filter on the 
scale enclosing mass M . Universality of the halo mass 
function has been tested in numerous simulations, with 
results generally confirming the assumption even if the 
specific functional forms for f{v) may differ from one 
another. 

The significance of a halo of mass M depends on 
the background density field so one can compute 
dn/d5i{x) and insert it into Equation (jH) [5l| . 



(16) 



(This is > 1 for massive haloes since the last derivative 
is negative in this case.) 

The derivative dhm/dlnag appearing in Equation 
(|14p can be obtained under the same universality assump- 
tion. In fact, the calculation is simpler. The definition 
of the significance implies v oc cTg"^, so that dlnv/dlnM 
does not depend on erg at fixed M. Therefore n oc i'/(i') 
and 

51nn ^ din. dH.fi.)] ^ _d 
am as am erg omi' dv 

Thus by comparison to Equation (|16p . we find: 



A6(M,A:) = 3/nl(6-1)<5c 



k^T{k)D{z) 



^] (18) 



This is equivalent to the previously derived expressions 
[47I [5^ . However, it is more general, because it is in- 
dependent of the form of f{v). It is therefore valid for 
the Press-Schechter mass function as well as for the more 
accurate Sheth-Tormen function. It is also valid for any 
object that obeys a halo occupation distribution (HOD) 
that depends only on the halo mass, {N){M), since in 
this case both b and Ab are linearly averages of their 
values for individual masses: 



/ b{M)n{M){N){M) dM 



J n{M){N){M)AM 



(19) 



and similarly for A6. 



The above statements apply to biasing of objects whose 
HOD depends only on the mass of the halo. However this 
may not be true for the quasars; in particular there are 
many lines of evidence that su gges ts that quasar activity 
is triggered by recent mergers |53l Is^. Issj . Therefore we 
should consider the standard bias b and large-scale bias 
A6 for recent mergers, which is in general not the same 
as the bias of all haloes of the final mass [56', '57| . 

This section considers the simplified case in which 
quasars are triggered by a merger between a halo of mass 
Ml and one of mass M2 , after which the quasar lives for a 
time tq <C in the new halo of mass Mq = Mi -f iV/2. 
This requires us to understand the dependence of the 
number of recent mergers on amplitude, which we will 
again express as erg- Unlike the case of the mass func- 
tion there are no accurate fitting formulae to the merger 
rate that have been tested against A^-body simulations 
for a variety of cosmologies. Therefore we will take two 
approaches here. The first will be to consider the recent 
merger probabilities from the extended Press-Schechter 
(ePS) formalism. With ePS, we will find that for haloes 
of a given mass the probability of being a recent merger 
is proportional to (J^^. In this picture, the bias of the 
quasars in this case is the same as the halo bias 6(Mo), 
but the /NL-induced bias A6 is less for recent mergers 
than for all haloes of mass Mq. However there is no rig- 
orous error bound on ePS calculations, so it is desirable 
to have an independent way to get the dependence of 
merger histories on We therefore consider a second 
method of getting the recent merger probability during 
the matter-dominated era by using redshift scaling rela- 
tions from A^-body simulation results. The latter method 
confirms the ePS relation. 



1. Extended Press-Schechter calculation 

We will work in the ePS formalism in which the merger 
history seen by a given dark matter particle is con- 
trolled by the linear density field 5{M) measured to- 
day, spherically smoothed on a mass scale M in La- 
grangian space [5^. At time t a particle is inside a 
halo of mass > M if 6{M') > uj{t) for any M' > M, 
where Lo{t) = 5cD{tQ)/ D{t) is the ratio of the threshold 
over-density for collapse Sc to the growth function D{t). 
In this picture it is convenient to replace the smooth- 
ing scale M with the variance of the density field on that 
scale, S{M) = {5{Mf). The smoothed density field 6{S) 
then follows a random walk as a function of S: this ran- 
dom walk is usually assumed to be Markovian because 
(i) each Fourier mode is independent for Gaussian initial 
conditions, and (ii) one neglects the difference between 
smoothing with a top-hat in Fourier space (in which each 
mode is independent) and the physically motivated top- 
hat in real space. 

In this formalism we would like the probability that a 
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halo of mass Mq at time t was actually of mass Mi at an 
earlier time t—tq and experienced a merger with a halo of 
mass A/2 — Mq — Ml. As argued by Lacey & Cole [ssl ]. 
the probability for a particular dark matter particle in 
this halo to have been in an object of mass < Mi at time 
t — tq is the probability that the trajectory 6{S) does not 
exceed uj{t — tq) between So = S{Mq) and Si ~ S{Mi). 
This evaluates to 

Ppa.ticie« Ml) = erfc (20) 



^2(^1 - ^0) 



so the differential probability is: 



Pparticlc(Mi)dMi 



1 Loit ~ tq) ~ uj{t) 



2^ (5i-5o)3/2 

[ioit - tq) - Luit)]^ 



X exp 



2(^1 - So) 



dSi 



dMi 



dMi 
tq\uj\ 



2^ {Si - 5*0)3/2 



dSi 



dMi 



dMi, 
(21) 



where in the last line we have assumed that tq is short 
so that one can do a Taylor expansion to lowest order in 
tq and recalled that a; < 0. The differential probability 
that the halo of mass Mo is a recent merger is simply this 
divided by the fraction of the particles in the mass Mi 
progenitor, 



P(Mi|Mo)dAfi 



1 



Mo 



dS-i 



dMi. 



(5i - 5o)3/2 Ml dMi 

(22) 

This formula was first derived by Lacey & Cole [58l |. 
albeit in a slightly different form [they computed 
P(Mo|Mi)]. It has two well-known deficiencies. One 
is that it is not symmetric under exchange of Mi and 
M2 , especially for extreme mass ratios [59,, . Another 
is that it does not contain merger bias in the Gaussian 
case, i.e. the bias of mergers is simply 6(Mo) [H^. This 
is because of the assumption that the trajectory 6{S) is 
a Markovian random walk, which is not quite correct. 
For example, the analytic explanations for merger bias 
of high-mass haloes 6l| are based on non-Markovian be- 
haviour due to the fact that the physically meaningful 
smoothing in real space does not correspond to a sharp 
cutoff at some fcmax- The corresponding merger history 
bias is based on the correlation coefficient 7 between 
6{S) and dS{S)/dS; one would have 7 = if one used 
the Fourier-space rather than real-space top-hat filter. 
This subtlety is however not required to understand why 
the merger bias in /nl cosmologies is significant on large 
scales. 

For our application we would also need to integrate 
over the range of masses Mi that define a major merger, 
but since the result does not actually depend on this we 
will not explicitly write it. In order to apply Equation 



p4)) to recent mergers we need only understand how the 
number density of recent mergers varies with ag. Since 
the number density of recent mergers is the product of the 
number density of haloes of mass Mo and the probability 
of them being recent mergers, we may write 



9 In n 



merger 



d\nn{Mo) ^ d\nP{Mi\Mo) ^^3) 
9 In (78 d In erg 



d In (jg 



From Equations (|T6l) and p7|) . the first term evaluates 
to 6c{b — 1). The second term contains the merger tree- 
dependent contribution to the large-scale bias, and can be 
evaluated from Equation (P^ . If one varies cg, the mass 
variances all scale as S{M) oc cr|, and hence P(Mi |Mo) oc 
cr^^. Thus we conclude that the second term is equal to 



-1, so 



9 In n-n 



d In (Tg 



(24) 



This is identical to the extra large-scale bias for haloes 
of fixed mass, except that we have a factor of 6 — 1 — (5^^ 
instead of 6 — 1, so the factor of 6 — 1 is replaced by 6 — 1.6. 

This formula is derived from ePS formalism and so 
it would seem to be on a somewhat less certain foot- 
ing, since the analytic formulas for merger rates have not 
been tested in A^-body simulations as extensively as those 
for the halo mass function. However, as we show in the 
following subsection, we do have some guidance from nu- 
merical simulations suggesting the scaling derived here is 
correct. 



2. Scaling from N-body simulations 

The ePS formalism predicts that the probability of a 
halo of mass Mo being a recent merger, P(Mi|Mo)dMi, 
is proportional to cr^^. While the qualitative result that 
massive haloes are more likely to be recent mergers in 
low-ag than high-ag cosmologies is supported by A^-body 
simulations [63 . [63| , the quantitative validity of the — 1 
exponent does not appear to be well-tested. Neverthe- 
less, in the matter-dominated era it is possible to deter- 
mine the exponent from the rcdshift dependence of the 
merger rate. 

The key is that there is no preferred timescale in 
the Einstein-de Sitter cosmology; scale factor and linear 
growth factor are both proportional to ^^/s and therefore 
the rescaling of initial amplitude is mathematically iden- 
tical to rescaling of time. Hence two A^-body simulations 
whose initial conditions differ only by the normalization 
of the primordial perturbations will evolve through ex- 
actly the same sequence of halo formation and mergers, 
except that the scale factor of each merger is re-scaled 

according to Omcrger oc CT^^. 

Note that for ACDM cosmologies this correspondence 
between scaling time and scaling normalisation breaks as 
the amplitude of fluctuations at the onset of cosmic ac- 
celeration will be different for different initial amplitudes. 
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Therefore results of this rescahng do not apply to the low- 
est redshifts, where the dark energy becomes important, 
but since in our analysis the only sample where recent 
mergers may be relevant is the quasar sample, which has 
a redshift distribution peaking at z ~ 1.7, this is a minor 
deficiency. 

We want to test this relation from the merger his- 
tory statistics [iJ] in the Millenium Simulation [6^. 
The recent merger probability, which we have denoted 
P(M-] |Mn)dMi , is related to the merger rate B/n de- 
fined by Ref. [6| by: 



P(Mi|Mo)dAfi =tQ 



n(Mo) 



(25) 



where ^ > 1 is the mass ratio of the progenitors. Here 
B /n is the merger rate per final halo of mass Mq per unit 
redshift per unit ^. The derivative with respect to erg is 
straightforward to express: 



(91nP(A/i|A/o) _ d\n{B/n) 
(9 In (Tg d In ag 



(26) 



where the partial derivatives are all at constant z. 

In an Einstein-de Sitter universe, the rescaling of the 
amplitude 

as ^ (1 + e)(T8 (27) 
is equivalent to rescaling of the scale factor (or redshift): 

1 + (28) 



1 + e 



Equating the recent merger probabilities in these two 
cases gives 



S , , , B fl + z 

— (1 + z, (78(1 + ejjdz = — [ — — , as 
n n \1 + e 



1 



(29) 



[The denominator 1-l-e on the right-hand side comes from 
rescaling of the redshift interval, dz — > dz/(l-|-e).] Taking 
the logarithm of both sides, and then differentiating with 
respect to e gives 



(91n(B/n) 
9 In as 



-1 



d\n{B/n) 
ain(l + z)' 



(30) 



This means that Equation (|26p can be evaluated provided 
the power-law exponent relating B/rt to 1 -I- z is known. 

The ePS prediction is B/n cx (1 -I- z)°, i.e. constant 
[S^ l- Inspection of Figure 8 of Ref. 164| shows that for 
a wide range of halo masses (> 2 x IQi^^Mq) and pro- 
genitor mass ratios (100:1 through 3:1), the exponent is 
indeed close to during the matter-dominated era z > 1, 
although in some cases (galaxy mass haloes, 3:1 mergers) 
the actual scaling is closer to B/n oc (l-|-z)°-^. These re- 
sults suggest the scaling 91nP(Afi|Mo)/91ncr8 is in the 
range of —1 (the ePS prediction) to —1.1. If we plug this 



into Equation 



then one derives (91nn„icrgor/i91n(T8 



These results provide an independent calculation of 
d\nP{Mi\Mo)/dhias that is on a more solid footing 
than cPS. The agreement of the logarithmic derivatives 
at the ~ 10% level is remarkable, especially given that 
ePS does not do so well at predicting the absolute merger 
rate. 



3. Summary 

We can write a generalised expression for the /nl in- 
duced scale dependent bias as 



A6(M,fc) = 3/nl(&-p)<5c 



nr. 



Ho 



k^T{k)D{z) \ c 



(31) 



equal to Sc{b - 1.6) (for -1) or 5 c{b - 1.65) (for -1.1) 



where 1 < p < 1.6, i.e. p — 1 for objects populating a 
fair sample of all the halos in a given mass range and 
p = l + 5~^ ~ 1.6 for objects that populate only recently 
merged halos. Below we discuss plausible values of p for 
the data samples used in this paper. 

To summarize, in non-Gaussian cosmologies, there are 
two types of merger bias: both b and A6 can depend 
on the merger history of a halo as well as its final mass. 
The ePS prediction for recent mergers is that for the 
bias b = b(Mo), i.e. there is no dependence on merger 
history; but that recent mergers with final mass Mq have 
a smaller Ab than one would find considering all haloes 
of mass Mq. Under the specific assumptions of ePS, if 
one makes the extreme assumption that all quasars are 
the result of recent halo mergers, the correction can be 
implemented by replacing & — 1 in Equation (jlSp with 
6-1.6. 

The reliability of the ePS result can only be evaluated 
by comparison to A^-body simulations. In the matter 
dominated era, in the range of masses and progenitor 
mass ratios covered by Ref. [gj], ePS appears to be a 
good description for the merger bias of Ab. 

Since in practice one estimates A6 from the observed 
clustering rather than from the unobserved halo mass 
Afo, any assembly bias effects in b (66j are also impor- 
tant to our analysis. This subject has received much 
attention recently, with the general result being that for 
high-mass haloes {M ^ A^*), those haloes that exhibit 
substructure, have lower concentration, or are younger 
have a slightly higher bias than the mean b{M). For 
example, in (67j it was found that the lowest quartile 
of haloes in concentration has bias ~ 10-20% higher 
than the mean b{M). Ref. [i^ found that the lowest- 
concentration quintile was 10% more biased than their 
highest-concentration quintile, and that this dependence 
was even weaker if one spli t on formation redshift instead 
of concentration. Ref. [69| found almost no dependence 
on formation redshift in the relevant range 8c/ a > 2, but 
their lowest quintile of concentration is 25% more bi- 
ased than the mean b{M) and evenj^arger effects are seen 
if one splits by substructure. Ref. [57] found that the bias 
for recent major mergers was enhanced by ~ 5% relative 
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to b{M). It is clear that the strength of this effect de- 
pends strongly on the second parameter used, but in the 
case of the definitions related to the mergers, the devia- 
tion in 6(M) is significantly less than the corrections to 
the Ab, which replaces 6—1 with 6 — 1.6. We will there- 
fore assume that the theoretical uncertainties are fully 
absorbed by the expression in Equation (j3ip . 

Finally, while there is ample evidence that quasar ac- 
tivity is often triggered by mergers, it is probably not 
the case that all quasars live in recently merged halos. 
Therefore the true value of p for quasar population lies 
between 1 and 1.6, since the true population of host ha- 
los lies somewhere between randomly selected halos and 
recently merged halos. Therefore, our limits with p ~ 1.6 
should be viewed as a most conservative reasonable op- 
tion. 



III. METHOD AND DATA 

After reviewing and extending the theoretical formal- 
ism we turn to its application to the real data. We would 
like to use Equation (|3ip to put constraints on the value 
of the /nl parameter. Since the effect is significant only 
on very large scales we need to use the tracers of large 
scale structure at the largest scales available. In addi- 
tion, the effect scales as b — p, where p is typically 1 but 
can be as large as 1.6 in special cases, hence we need very 
biased tracers of large scale structure to measure the ef- 
fect. We discuss below our choice of observational data. 
Finally, the effect changes the power on large scales and 
in principle this can also be achieved with a change in 
the primordial power spectrum, although this degener- 
acy exists only in the presence of one tracer: with two 
tracers with different biases one can separate /nl from 
the changes in the initial power spectrum. Here we as- 
sume that the basic model is one predicted by the sim- 
plest models of structure formation and we do not allow 
for sudden changes in the power spectrum beyond what 
is allowed by the standard models, which assume the 
power spectrum slope Ug to be constant. We use Markov 
Chain Monte Carlo method to sample the available pa- 
rameter space using a modified version of the popular 
public package cosmomcjTOj. In addition to /nl, we fit 
for the standard parameters of the minimal concordance 
cosmological model: ujh = flhh'^, i^cdm = f^CDM^^, 0, t, 
Us and logj4, where 9 is the ratio of the sound horizon to 
the angular diameter distance at decoupling (acting as a 
proxy for Hubble's constant), r is teh optical depth and 
A is the primordial amplitude of the power spectrum. 
All priors are wide enough so that they do not cut the 
posterior at any plane in the parameter space. 

We always use standard cosmological data as our base- 
line model. These include the WMAP 5-year power- 
spectra [tiI, [T^I and additional smaller-scale experiments 
(VSA, CBI, ACBAR) [zl, H Izl, as well as the super- 
novae measurements of luminosity distance from the Su- 
pernova Legacy Survey (SNLS) For values of /nl 



under consideration in this paper these data sets are not 
directly sensitive to the /nl parameter. However, they 
are needed to constrain the basic cosmological model and 
thus the shape and normalisation of the matter power 
spectrum. The large scale structure data discussed bel- 
low are thus simultaneously able to fit for /nl and other 
auxiliary parameters. 

Most of our large-scale structure data is drawn from 
the Sloan Digital Sky Survey (SDSS)^ The SDSS drift- 
scans the sky in five bands (ugriz) [73] under photomet- 
ric conditions [tIjItI] using a 2.5-meter optical telescope 
[80| with 3 degree field of view camera [8l| located in New 
Mexico, USA [TSj. The photometric and astrometric cal- 
ibration of the SDSS and the quality assessment pipeline 
are described by Refs. [85, "ss", "si, 8§, [sl. Bright galax- 
ies jSTi], luminous red galaxies (LRGs) [88[, and quasar 
candidates [8§| are selected from the SDSS imaging data 
for spectroscopic follow-up ^9^. This paper uses imag- 
ing data through the summer of 2005, which was part of 
SDSS Data Release 6 fDR 6') iQll . and spectroscopic data 
through June 2004 (DR4) [ofT^ 

The requirement of large scales and highly biased trac- 
ers leads us to explore several different large scale data 
sets, in particular SDSS LRGs, both spectroscopic 
and photometric 193 samples, and photometric quasars 
(QSOs) from SDSS [9g|- In all these cases we use auto- 
correlation power spectrum, which is sensitive to /nl- 
These datasets can be assumed to be statistically inde- 
pendent. As explained below, when spectroscopic and 
photometric LRGs are analysed together, we take care 
exclude those photometric redshift bins that have sig- 
nificant overlap with spectroscopic sample. Our quasars 
have typical redshifts of z ~ 1.5-2 and only 5% overlap 
with LRGs due to photoz errors. Moreover, they are in 
the Poission noise limited regime, so overlap in volume is 
less relevant. This is discussed in more detail in [9^1 • In 
addition, we also use cross-correlation of all these sam- 
ples, as well as nearby galaxies from the 2-micron All-Sky 
Survey (2MASS) [sg and radio sources from the NRAO 
VLA Sky Survey (NVSS) [HI with CMB, as analyzed in 
Ho et al. [9^. Since this is a cross-correlation between 
the galaxies and matter (as traced by the ISW effect), 
the dependence is linear in /nl- 

A. Spectroscopic LRGs from SDSS 

We use the spectroscopic LRG power spectrum from 
Tegmark et al. [93| , based on a galaxy sample that covers 
4000 square degrees of sky over the redshift range 0.16 < 
z < 0.47. We include only bins with k < 0.2/i/Mpc. We 
model the observed data as 

Pobscrved(A:) = [b + Ab{k, /nl)]' ^lin (fc) ^^^^^^ , (32) 

where the last term describes small-scale non-linearities 
[131 , with Q treated here as a free parameter (bound be- 
tween zero and 40) and A = lAh~^ Mpc. For realistic 
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values of /nl and Q, the non-Gaussian bias Ab is present 
only at large scales and the Q-term is present only at 
small scales; there is no range of scales over which both 
are important. We explicitly confirmed that there is no 
correlation between Q and /nl present in our MCMC 
chains and we let the data to determine the two param- 
eters. 

As discussed above if the halos in which objects re- 
side have undergone recent mergers and thus depend on 
properties other than halo mass then the simple scaling 
with (b — 1) may not be valid. This is unlikely to be rel- 
evant for LRGs, which are old red galaxies sitting at the 
center of group and cluster sized halos. Moreover, num- 
ber density of LRGs is so high that it is reasonable to 
assume that almost every group sized halo contains one, 
since otherwise it is difficult to satisfy both the number 
density and high bias at the same time [98,]. For LRGs 
we thus do not expect there is a second variable in addi- 
tion to halo mass and halo occupation models find that 
populating all halos with mass above IO^^Mq/Zi with an 
LRG is consistent with all the available data [93 • Hence 
we will only use p = 1 in Equation ([3T|) . The same also 
holds for the photometric LRGs discussed below. 



B. Photometric LRGs from SDSS 

We use data from Padmanabhan et al. [o^, who pro- 
vide the LRG angular power spectrum measured in 8 
redshift slices (denoted 0-7) covering the range 0.2 < 
2^photo < 0.6 in slices of width Azphoto = 0.05. The power 
spectrum is based on 3500 square degrees of data. We 
use only data-points that correspond to fc < 0.1/i/Mpc 
at the mean redshift of each individual slice. We use 
the full Bessel integration to calculate the angular power 
spectrum on largest scales and account for the redshift- 
distortion power as described in Padmanabhan et al. [o^l : 



Ci = cr + cr + c, 



vv 



(33) 



where superscripts g and v denote galaxies over-density 
and velocity terms respectively. The bias and /3 depen- 
dence has been put back into the Bessel integral as it now 
depends on the value of fc. The three terms are given by 
the integrals: 



A\k)\W,[k)\\ 



r<99 


^4, J 


f dfc 

T 
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f dfc 
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where A^(fc) is the linear matter power spectrum today. 
The window functions are given by 



W[{k) 



D{0) dr 
D{r) dN 



and 



D{0) dr 



2£2 -f 2^ - 1 



(2^- 1)(2£-|- 1) 
{£+!){£ + 2) 
' {2e+l){2e + 3) 



{2£^1){2£- 
je-2{kr) 



3) 



dr, 



(35) 



(34) 



where dN/dr is the redshift distribution normalised to 
unity and written in terms of comoving distance r, and D 
is the growth function. The code automatically switches 
to the Limber approximation when this becomes accu- 
rate. Note that for the low multipoles, it is essential to 
include the redshift-space distortion even for a photomet- 
ric survey because a significant amount of power comes 
from Fourier modes that are not transverse to the line of 
sight. 

We use an independent bias parameter for each red- 
shift slice. In addition to the bias dependence, we also 
use the Equation (j32p to take into account non- linear cor- 
rections. While strictly speaking the value of Q should be 
different for each slice, we use a single free parameter Q 
for all slices. We have explicitly checked that non-linear 
corrections are negligible for fc < 0.1/i/Mpc and there- 
fore this is not a major issue. Since there is a strong 
overlap between the spectroscopic and photometric sam- 
ple for z < 0.45, we use only slices 5-7 when combining 
this data with spectroscopic sample. 



C. Photometric quasars from SDSS 

We use the power spectrum of the high redshift 
quasar photometric sample that has recently been con- 
structe d for ISW and CMB lensing cross-correlation stud- 
ies [il, llOOj l . The sample covers 5800 square degrees of 
sky and originally had two photometric redshift ranges, 
0.65 < zphoto < 1.45 ("QSOO") and 1.45 < ^photo < 
2.00 ("QSOl"). It consists of ultraviolet-excess (UVX; 
u — g < 1.0) point sources, classified photome trica lly as 
quasars [lOl( , and with photometric redshifts [lO^ ■ The 
classification and photometric redshifts were at the time 
of sample construction only available over a subset of the 
survey region; they were extended to the remaining re- 
gion using a nearest-neighbor algorithm in color space 
|95j . As described below, we only used the QSOl sam- 
ple as QSOO appears to suffer from systematic errors on 
large scales. 

The largest angular scales in the quasar data are sub- 
ject to at least three major sources of systematic error: 
stellar contamination, errors in the Galactic extinction 
maps, and calibration errors. All of these are potentially 
much worse than for the LRGs: some stars (e.g. M dwarf- 
white dwarf binaries) can masquerade as quasars, and we 
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are relying on the u band where extinction is most severe 
and the photometric cahbration is least well understood. 
These errors were discussed in the context of ISW stud- 
ies [95*1, but if one wishes to use the quasar auto-power 
spectrum on the largest angular scales the situation is 
more severe. 

We investigated this subject by computing the cross- 
power spectrum of each QSO sample with the SDSS 
18.0 < r < 18.5 star sample and with "red" stars (which 
satisfy the additional cut g — r > 1.4). The cross-power 
should be zero in the absence of systematics but it could 
be positive if there is stellar contamination in the QSO 
sample. Either positive or negative correlation could re- 
sult from photometric calibration errors which shift the 
quasar and stellar locus in nontrivial ways. The results 
of this correlation are shown in Figure [1] with error bars 
estimated from the usual harmonic space method, 



QSOO-star power spectrum 



QS01 -star power spectrum 



^(Cr + n,-i)C| 



1 



^min] /sky 



(36) 



where C'' is the quasar autopower spectrum, Ug is the 
number of quasars per steradian, and C'l'* is the star au- 
topower spectrum; the nj^ term is negligible. (Aside 
from boundary effects, this is the same error that one 
would obtain by correlating random realizations of the 
quasar field with the actual star field.) From the fig- 
ure, QSOl appears clean, but QSOO appears contami- 
nated: the first bin {2 < £ < 12) has a correlation of 
= -(2.9 ± 1.0) X 10-4 with the red stars. This is 
a — 2.9ct result and strongly suggests some type of sys- 
tematic in the QSOO signal on the largest angular scales. 
We are not sure of the source of this systematic, but the 
amount of power in the quasar map that is correlated 
with the red stars is 



2tt 



(corr) 



^(^+1) C, 
2-K 



qs 2 



2 X 10" 



(37) 



The variation in QSOO density that is correlated with the 
red stars is thus at the ^ 1.4% level. This is consistent 
with the excess auto-power in QSOO in the largest scale 
bin, which is at the level of [£{£ + l)/27r]C^''' - 3 x lO""^ 
(see figure 10 in 93|) and is comparable to what one 
might expect from the 1-2% calibration errors in SDSS 
[84| . although other explanations such as the extinction 
map are also possible. Because of this evidence for sys- 
tematics, we have not used the QSOO autopower spec- 
trum in our analysis; in what follows wc only use QSOl. 
QSOO may be added in a future analysis if our under- 
standing of the systematics improves. 

Because the redshift distribution is poorly known and 
needs to be determined internally from the data sample 
itself the data is analyzed in a two step procedure: 

1. Power spectrum points with £ = 30 . . . 200 are used 
to constrain the product {bdn/dz){z), where b is the 
linear bias at redshift z and dn / dz is the normalised 
radial window function. Although the /nl is taken 
into account at this step, its effect is sub-dominant. 
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FIG. 1: The correlation of each quasar sample with stars and 
with the red stars. 



2. We then assume that the functional form for b{z) 
is either 



b{z) (X.I + 



l + z 
2.5 



(38) 



as measured in |103j , or that its form is given by 

b{z) ozl/D{z), (39) 

as would be valid if the clustering amplitude is not 
changing with redshift. In both cases the constant 
of proportionality is determined from the normali- 
sation condition 



dn 
dz 



-dz 



bdn 1 
"d76(^)' 



-dz = 1. 



(40) 



3. Once both dn/dz and b{z) are known we calculate 
the theoretical angular power spectrum using the 
same code as for photometric LRGs, taking into 
account all available £ points. This is the theoreti- 
cal spectrum that is used to calculate the that 
goes into the MCMC procedure. In principle, all 
free parameters that determine the shape of bdn/dz 
should be varied through the MCMC, rather than 
being fixed at the best-fit point. However, in the 
limit of Gaussian likelihood, the two procedures are 
equivalent, while the latter offers significant speed 
advantages. 

4. We calculate using two different methods. Our 
standard method is to use all points and the full 
covariance matrix assuming a Gaussian likelihood: 



X 



(d-t) •C-i(d-t), 



(41) 



where d and t are data and theory vectors of 's, 
respectively. This is likely to be a good approxima- 
tion except on the largest scales, where small num- 
ber of modes leads to corrections that may affect 
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QSO autopower window functions 




FIG. 2: The window functions for the £ — 6 and 
of the quasar power spectrum. 



= 18 bins 



the outer limits (e.g. 3a). To test this we model 
the first QSO bin with an inverse-x^ distribution, 
but neglecting covariance of this bin with higher £ 
bins: 



= N 



In 



di + l/n\ ti + l/n 



te + l/n J di + l/n 



- 1 



Gaussian for other points, (42) 



D. Cross-correlation between galaxies and dark 
matter via Integrated Sach- Wolfe effect 

One can also look for the /nl using the cross- 
correlation between a tracer like galaxies or QSOs and 
dark matter. Since we do not have dark matter maps 
from large scales we can use cosmic microwave back- 
ground (CMB) maps as a proxy. If the gravitational 
potential is time dependent, as is the case in a universe 
dominated by dark energy or curvature then this leads 
to a signature in C MB, the so-called integrated Sach- 
Wolfe (ISW) effect |l05j . This signature can easily be 
related to the dark matter distribution, but part of the 
signal is not coming from ISW but from the primary 
CMB anisotropics at the last scattering surface. These 
act as a noise and lead to a large sampling variance on 
large scales and as a result the statistical power of this 
technique is weaken ed. A t the moment ISW is only de- 
tected at ^ 4(7 level jlQ6l | . Our procedure closely follows 
that of [9^ and is in many respects very similar to that 
of the Section lill CI We use all 9 samples present in Ho 
et al., although the discriminating power is mostly com- 
ing from the NVSS-CMB cross-correlation because NVSS 
sample is available over 27 361 deg^ area and the tracers, 
radio galaxies, are biased with 6^2. First, bdn/dz and 
b{z) are determined for each sample. Here we always use 
b{z) oc 1/D{z) for all samples except the quasar sample. 
Then we calculate the ISW Limber integral: 



where N ~ 21.5 is the effective number of modes 
contributing to the first bin and n is the number of 
quasars per steradian. (Thi s was called the "equal 
variance" case in Ref. (104| . and is appropriate if 
there is equal power in all modes; this is the case 
here to a first approximation since we find that the 
Poisson noise l/n dominates.) 

Since the power rises dramatically at low i in the /ml 
models, we have included the full window function in 
our calculation of the binned power spectrum, Cbin = 
WeCe. The window functions for the lowest two QSO 
bins are shown in Figure [21 and the algorithm for their 
computation is presented in Appendix [XI 

As shown in the Section [II C li if halos have undergone 
recent mergers then we would expect p 1.6 instead of 
p ~ 1 in Equation (|3ip . Recent mergers could be a plau- 
sible model for QSOs, whose activity could be triggered 
by a merger [HI, HI, |5l| . For QSOs used in this work we 
do not find that their number density places a significant 
constraint: at z ~ 1.8 the number density of halos with 
6 2.5 — 3 is one to two orders of magnitude higher than 
the measured number density, hence we can easily pick 
and choose the halos with a recent merger and still satisfy 
the combined number density and bi as c onstraint. This 
is in agreement with conclusions of [l03l |. which looked 
at a similar QSO sample from 2dF. Because of this un- 
certainty we therefore run another separate analysis for 
QSOs with p = 1.6 (QSO merger case). 



^ C3(£+ 1/2)2 



dz [b{z) + Ab{k{z),z)] 



dn d 
dz dz 



m 



(1 + ^) 



Diz)P{kiz)), (43) 



where k{z) — {£ + 1/2) /x. Again, we used the full win- 
dow function for the NVSS-CMB correlation to avoid an 
unnecessary bias and we assume p = 1 for all the samples. 



IV. RESULTS 

We begin by plotting data-points and theoretical pre- 
dictions for six of the datasets and values of /nl in Figure 
[31 This plot deserves some discussion. The easiest and 
most intuitive to understand is the case of spectroscopic 
LRGs. We note that the inclusion of the /nl param- 
eter modifies the behaviour of the power spectrum on 
the largest measured scales. Naively, one would expect 
LRGs to be very competitive at constraining /nl- In 
practice, however, /nl is degenerate with matter den- 
sity, which also affects the shape of the power spectrum 
and hence the constraints are somewhat weaker. Since 
the effect of /nl raises very strongly, just a couple of 
points on largest scales might break this degeneracy, im- 
proving constraints by a significant factor. The effect 
in the photometric LRG samples is similar, although we 
are now looking the angular space, where the dependence 
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FIG. 3: This figure shows 6 datasets that are most relevant for our constraints on the value of /nl- In the left column show 
the NVSSxCMB Integrate Sach Wolfe Cross correlation, the QSOl power spectrum, the spectroscopic LRG power spectrum, 
while the right column shows the last three slices of the photometric LRG sample. The lines show the best fit /nl = model 
(black, solid) and two non-Gaussian models: /nl ~ 100 (blue, dotted), /nl ~ —100 (red, dashed). The ISW panel additionally 
shows the /nl = 800 model as green, dot-dashed line. While changing /nl, other cosmological parameters were kept fixed. See 
text for further discussion. 
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Phot LRG 
Phot LRG (0-4) 
Spec LRG 
ISW 
QSO 
QSO (b=1/D) 

QSO alt 
QSO merger 
Combined 
Comb, merger 
Comb.+W5 bi. 
Merg.+WS bi. 
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400 



FIG. 4: This figure shows the median value (red points) and 1,2 and 3-sigma Umits on /nl obtained from diflerent probes 
(vertical lines). The data set used are, from top to bottom: Photometric LRGs, Photometric LRGs with only slices 0-4 used. 
Spectroscopic LRGs, Integrated Sach Wolfe effect, photometric QSO, photometric QSOs using b[z) oc 1/D(z) biasing scheme 
(see Section Till C[) . photometric QSOs using alternative calculation scheme (see Section Fill Cfl . using a scale dependent bias 
formula appropriate for recently merged halos (Section III C|) . Combined sample. Combined sample using a scale dependent bias 
formula appropriate for recently merged halos (for QSO), the last two resoluts to which a statistically independent WMAP 5 
bispectrum /nl constraint was added. See text for discussion. 



has been smeared out. The QSO plot again shows simi- 
lar behaviour, with two caveats. First, the changes in the 
predicted power spectrum on smaU scales are a result of 
the fact that hdn/dz is perturbed with changing /nl, al- 
though this is a minor effect. Second, the increase in the 
power at smallest £ for negative /nl is due to the fact that 
for sufficiently negative /nl (or sufficiently large scales), 
5h < —2b and hence power spectrum rises again above 
what is expected in the Gaussian case. The more unex- 
pected is the NVSS-CMB cross-correlation. Naively, one 
would expect that the first point of that plot will produce 
a very strong /nl "detection" . However, the CMB cross- 
correlation signal is only linearly dependent on /nl, while 
cross-correlations of NVSS with other tracers of struc- 
ture are quadratically dependent on /nl. Large values 



of /nl produce anomalously large power in the angular 
power spectrum if bdn/dz has significant contributution 
at high-z tail, which probes large scales. Therefore, the 
6dn/dz fitting procedure skews the distribution towards 
lower redshifts, leading to a lower bias overall. At very 
large values, e.g. /nl — 800, this effect is so severe that 
the b oc 1/D{z) scaling forces 6 < 1 at the low-redshift 
end. This implies Ab < 0, and the large-scale ISW sig- 
nal actually goes negative (see top-left panel of Fig. [3]). 
Therefore, the ISW is surprisingly bad at discriminating 
/nl and we were unable to fit the first NVSS ISW data 
point with a positive /nl. This behavior is however only 
of academic interest because the other data sets strongly 
rule out these extreme values of /nl. 

We ran a series of MCMC chains with base cosmolog- 
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Data set 



Photometric LRG 
Photometric LRG (0-4) 
Spectroscopic LRG 
ISW 
QSO 
QSO (b=l/D) 
QSO ahernative ^ 
QSO merger 
Combined 
Comb, merger 
Combined + WMAP5 bispectrum /nl 
Combined merger + WMAP5 bispectrum /nl 



/nl strains /nl in both directions. The second and third 

point have some excess power relative to the best fit 
/nl = model and so give rise to a slightly positive value 
ii5+2i5+3ooof /nl in the final fits. This is however not a statistically 
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TABLE I: Marginalised constraints on /nl with 
% c.L), 2a (95% c.l.) and 3cr (99.7%) errors. 



mean, Icr 
Ordering 

exactly matches that of the Figure |4] See text for discussion. 



ical data and one of the four data sets considered above, 
as well as a run in which all data were combined (with 
the exception of slices 0-4 of photometric LRGs as de- 
scribed above) . The results are summarized in Table |T] 
and visualized in the Figure HI 

We note several interesting observations. When only 
slices 0-4 of the photometric LRG sample are used, the 
results are weaker than those of the spectroscopic LRGs. 
The two trace comparable volume, but photometric sam- 
ple does not use radial mode information and as a re- 
sult its errors are expectedly larger. On the other hand, 
the overall photometric sample performs somewhat bet- 
ter than the spectroscopic LRG sample, due to its larger 
volume: it traces LRGs up to 2; ~ 0.6 as opposed to 
z ~ 0.45 for spectroscopic sample. 

We find that the ISW is constraining the /nl param- 
eter rather weakly. This is somewhat disappointing, but 
not surprising, since the cross-correlation between LSS 
and CMB is weak and has only been detected at a few 
sigma overall. In addition, as mentioned above, /nl en- 
ters only linearly (rather than quadratically) in the ISW 
expressions and is strongly degenerate with determina- 
tion of hAnj^z. Given that ISW SjN can only be im- 
proved by another factor of ~ 2 at most even with perfect 
data we do not expect that it will ever provide competi- 
tive constraints on /nl- 

We find that the quasar power spectra give the 
strongest constraints on /nl- This is due to their large 
volume and high bias. The lowest I point in this data 
set at £ = 6 is a non-detection and therefore highly con- 



significant deviation from /nl = 0. 

We also test the robustness of /nl constraints from 
^g^^_i_^gg^gj2 X^oR quasar sample power spectrum by performing the 
i28following tests. First we replace the form of evolution of 
bias from that in Eq. (f38|) to that in Eq. (|39|) . Second we 
replaced the details of the likelihood shape for the first 
points, we replace Eq. (|^T|) with Eq. (|^^ . In both cases, 
the effect on the constraints was minimal, as shown in 
the Figure Hand the Table H 

The large-scale quasar power spectrum could also be 
affected by spurious power (e.g. calibration fluctuations 
or errors in the Galactic extinction map) or 1-halo shot 
noise, either of which would dominate at large scales over 
the conventional P(fc) oc k autocorrelation. At the level 
of the current data these are not affecting our results: 
uncorrelated power adds to the power spectrum Q , and 
if present it would only would only tighten the upper 
limits on the power spectrum given by our £ = 6 and 18 
points. The 1-halo shot noise would contribute an added 
constant to Ci since z ~ 1.7 haloes are unresolved across 
our entire range of l\ the observed power at £ ~ 250 
constrains the 1-halo noise to be much less than the error 
bars at £ = 6, 18. Nevertheless, if we had detected excess 
power in the quasars at the largest scales, a much more 
detailed analysis would have been necessary to show that 
it was in fact due to /nl and not to systematics. 

If quasars are triggered by mergers, then they do 
not reside in randomly picked halos and the use of the 
Eq. (fT8|) may not be appropriate. As discussed in Sec- 
tion [ircl one can replace the (6 — 1) factor in Eq. (fT8|) 
with a modified factor (h — 1.6). We quote the results 
in this case in the Table U as "QSO merger." The error 
bars have increased by about 40%, which is somewhat 
less than naively expected from a population with an av- 
erage bias of around 2.5. This is due to various feedback 
related to hAnj^z fitting. Moreover, this is hkely to be 
an extreme case since it assumes that all quasars live in 
recently merged halos with short lifetimes, so we would 
expect that the true answer is somewhere in between the 
two results. On the other hand, it is also unclear how 
accurate extended Press- Schechter formalism is for this 
application, so there is some uncertainty associated with 
this procedure. Recent progress in understanding quasar 
formation and its relation to the underlying halo popu- 
lation [54] gives us hope that this can be solved in the 
future. 

We also quote results for the combined analysis. In this 
case, the error on /nl shrinks somewhat more than one 
would expect assuming that error on /nl measurement 
from each individual dataset is independent. This is be- 
cause other parameters, notably matter density and am- 
plitude of fluctuations get better constrained when data 
are combined. Furthermore, we note that when quasars 
are assumed to be recent mergers the errors expand by 
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an expected amount at the Scr level, but hardly at a 2(t 
level. We have carefully investigated this anomaly and it 
seems to be due a peculiar shape of the likelihood surface 
and subtle interplay of degeneracies between /nl, matter 
density and amplitude of fluctuations. 



V. DISCUSSION 

The topic of this paper is signature of primordial non- 
Gaussianity of local type (the so called /nl model) in the 
large scale structure of the universe. Specifically, it was 
recently shown that this type of non-Gaussianity gives 
rise to the scale dependence of the highly biased trac- 
ers [13 . We extend this analysis by presenting a new 
derivation of the effect that elucidates the underlying as- 
sumptions and shows that in its simplest form it is based 
on the universality of the halo mass function only. Our 
derivation also allows for possible extensions of the sim- 
plest model, in which tracers of large scale structure de- 
pend on properties other than the halo mass. One that 
we explore in some detail is the effect of recent merg- 
ers. By using extended Press-Schechter model we show 
that in this case the effect has the same functional form, 
but the predicted scaling with /nl for a given bias has a 
smaller amplitude than in the simplest model. 

Second, we apply these results to constrain the value 
of /nl from the clustering of highly biased tracers of 
large scale structure at largest scales. In our analysis we 
find that the best tracers are highly biased photometric 
quasars from SDSS at redshifts between 1.5 and 2, fol- 
lowed by photometric and spectroscopic LRGs at redshift 
around 0.5. Our final limits at 95% (99.7%) confidence 
are 



29 (-65) < /nl < +70 (-f 93) 



(44) 



if we assume halos in which tracers reside are a fair sam- 
ple of all halos of a given bias. If we assume instead that 
QSOs are triggered by recent mergers and have short life- 
times then we find 



-31 (-96) < /nl < +70 (+96), 



(45) 



a somewhat weaker, but still competitive constraint. In 
both cases, we find no evidence for non-zero /nl- These 
results shows that existing data can already put very 
strong limits on the value of /nl, which are competi- 
tive with the best constraints from WMAP 5 year anal- 
ysis of CMB bispectrumQ, which is statistically inde- 
pendent of the method used in this paper. These give 
—9 < /nl < 111 at 95% confidence. Assuming that 
WMAP 5 constraint on /nl to independent of and well 
described by a Gaussian likelihood /nl = 51 ± 310, we 
get the following combined constraint: 



(-21) < /nl < +69 (+88). 



(46) 



If we assume quasars to live in recently merged halos, 
we get essentially the same result with the upper limit 
relaxed to 89. 




0.90 0.92 0.94 0.96 0.98 1.00 1.02 

Us 



FIG. 5: This figure shows 1 and 2a contours on the jIs-Jne 
plane for our best combined data set with additional WMAP 
5 year bispectrum constraint (assumed to be independent of 
Us). Red lines are predictions from the ekpyrotic models and 
correspond to values of fixed 7 and varying e. Different line 
correspond to different values of 7, which varies between 7 = 
— 1 (flat, constant, negative /nl) to 7 = —0.2 in steps of 0.1. 
The dashed green line corresponds to th e theoretically favored 
value of 7 = — 1/\/3 according to [107|]. 



In this combined result, /nl = is at just around 2a, 
which taken at a face value suggests than the evidence 
for a significant non-Gaussianity found in 3 year WMAP 
data by [? ] may have been a statistical fluctuation 
rather than evidence of a real signal. 

The results are already sufficiently strong to con- 
strain the ekpyrotic models of generating the primor- 
dial structure: these generically predict much higher 
non-Gaussianity than infiationary models 32]. Following 
Ref. [s^l, we parametrize the ekpyrotic model in terms 
of 7 = 02 /'/'I during the phase of creation of entropic 
perturbations for minimally coupled fields (/>i responsible 
for ekpyrosis and ekpyrotic parameter Cck S> 1 , the ekpy- 
rotic equivalent of the slow-roll parameter e. In Figure 
[5] we show our "Combined" constraints on the Hs - /nl 
plane together with approximate theoretical predictions 



/i 



NL 



2 1 logcck 

Cek 60 log 60' 

Eck — ot). 

7 



(47) 
(48) 



We note that large parts of parameter space are strongly 
constrai ned. For the theoretically preferred value of 7 = 
— 1/V3 (l07l |. eek is constrained to be between 100 and 
10000. Higher values of eck require considerably lower 
values of 7 and 7 > —0.3 is disfavored at 2a. 

Our results are very promising and with this first anal- 
ysis we already obtain constraints comparable to the best 
previous constraints. However, we should be cautious 
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and emphasize that this is only the first apphcation of 
this new method to the data and there are several is- 
sues that require further investigation. While analyti- 
cally the method is well motivated and can be derived 
on very general grounds, as shown in Section III C 1[ the 
method needs to be verified further in iV-body simula- 
tions using large scale tracers comparable to these used 
in our analysis. Equation (jlSp has been calibrated with 
iV-body simulations using matter-halo cross-correlation 
[13] ■ Auto-correlation analysis of biased halos, which is 
the basis for the strongest constraints derived here, is 
typically noisier and has not been verified at the same 
level of accuracy, although there is no obvious reason 
why it should give any different results. Still, it would be 
useful to have larger simulations where the scale depen- 
dent bias could be extracted with high significance from 
auto-correlation analysis. In addition, it would be useful 
to verify the scaling relations in simulations on samples 
defined as closely as possible to the real data should be 
selected, in our case by choosing halos with mean bias of 
2 at 2 = 0.5 for the LRG and 6 = 2.7 at z = 1.7 for QSO 
samples. 

A second uncertainty has to do with the halo bias de- 
pendence on parameters other than halo mass. We have 
shown that for QSOs we need to allow for the possibility 
that they are triggered by recent mergers and we have 
presented extended Press Schechter predictions for the 
amplitude of the effect in this case, which can affect the 
limits. To some extent these predictions have been ver- 
ified using simulations (63 |. but we do not have a re- 
liable model of populating QSOs inside halos to predict 
which of the limits is more appropriate for our QSO sam- 
ple. Moreover, just as the overall halo bias has recently 
been shown to depend on variables other than the halo 
mass [6^ , it is possible that the large scale /nl correction 
also depends on variables other than the mass and merg- 
ing history: while we have shown that extended Press- 
Schechter already predicts the dependence on the merg- 
ing history, other dependences less amenable to analytic 
calculations may also exist. It is clear that these issues 
deserve more attention and we plan to investigate them 
in more detail in the future using large scale simulations 
combined with realistic quasar formation models. 

On the observational side, there are many possible ex- 
tensions of our analysis that can be pursued. In this 
paper we have pursued mostly quasars and luminous red 
galaxies (LRGs) , two well studied and highly biased trac- 
ers of large scale structure. On quasar side, we have 
only analyzed photometric QSO samples split by redshift, 
but we should be able to obtain better constraints if we 
also use the luminosity informatio n, sp ecially if brighter 
quasars are more strongly biased [l03j |. Moreover, it is 
worthwhile to apply this analysis a lso t o > 3 spec- 
troscopic quasar sample from SDSS jl08l |. which is very 
highly biased. Its modelling appears to require almost 
every massive halo to host a quasar |l09l |. which would 
reduce the uncertainties related to secondary parameter 
halo bias dependence. An order of magnitude estimate 



shows that small values of /nl do not change mass func- 
tion enough to affect this deduction. On the LRG side, 
the most obvious extension of our work would be to pur- 
sue the luminosity dependent clustering analysis of the 
spectroscopic sample. It is well known that LRG cluster- 
ing amplitude is luminosity dependent jllOl | and so select- 
ing only brighter LRGs would lead to a higher bias sam- 
ple and could improve our limits, but ideally one would 
want to perform the analysis with optimal weighting to 
minimize the large scale errors as in [93| . Similar type of 
luminosity dependent analysis could also be done on the 
LRG photometric sample used here (o^ . 

Finally, with the future data sets from several planned 
or ongoing surveys, especially those related to the bary- 
onic oscillations most of which use highly biased tracers 
of large scale structure, a further increase in sensitivity 
should be possible [43|. The ultimate application of the 
large-scale clustering method would involve oversampling 
the 3-dimensional density field in several samples with a 
range of biases, so that excess clustering due to /nl can 
be cleanly separated from contamination due to errors 
in calibration or extinction correction, which are major 
challenges as one probes below the 1% level. The ulti- 
mate sensitivity of the method will likely depend on our 
ability to isolate these systematics and these should be 
the subject of future work. 

Overall, the remarkably tight constraints obtained 
from this first analysis on the real data, while still sub- 
ject to certain assumptions, is a cause of optimism for 
the future and we expect that the large scale clustering 
of highly biased tracers will emerge as one of the best 
methods to search for non-Gaussianity in initial condi- 
tions of our universe. 
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The parameters are then estimated according to 

Pi = {F'^^)tj{<l3 - ('7j>noisc), (A4) 

where (gj) noise is the expectation value of qj for Poisson 
noise. This is equal to 

(9j)noisc - I Tr (C-iC,C-iN) , (A5) 

where N is the Poission noise matrix (i.e. a diagonal 
matrix with entries equal to the reciprocal of the mean 
number of galaxies per pixel) . This quantity can be com- 
puted using the same machinery as used to compute F. 

The actual expectation values of the binned power 
spectra pi for a general power spectrum [i.e. not nec- 
essarily Equation (jAip ] are given by 

{Pi) ^ {F^^)i3{{q]) total - (gj)noisc); (A6) 

since signal and noise are uncorrelated this reduces to: 



APPENDIX A: WINDOW FUNCTIONS 

This appendix describes the computation of the win- 
dow functions. We begin with a brief revisit of the prin- 
ciples underlying the window function (see the references 
for more details) and then describe our computational 
method. 



1. Principles 



The power spectra in this paper were cor nput ed using 
the methodology of Padmanabhan et al. Ill ill, i mple- 



mented on the sphere as described in Refs. [94, 

urn- The 

power spectrum is estimated from a vector x of length 
TVpix containing the galaxy overdensities in each of the 
A^pix pixels. Following the notation of Ref. [93|, we esti- 
mate the power spectrum in bins. 



(Al) 



where is 1 if multipole I is in the ith bin, and other- 
wise, and Pi are the parameters to be estimated. We can 
then define the template matrices, C^, which are the par- 
tial derivatives of the covariance matrix of x with respect 
to Pi . The quadratic estimators are then defined: 



= ix^C-iQC-ix, 



(A2) 



where C is an estimate of the covariance matrix used to 
weight the data (our choice is described in Ho et al. [95|). 
We also build a Fisher matrix. 



1 



Tr (C-iQC-^C, 



(A3) 



(ft) = (^^"i).,(g,).ignai = ^iF-%TT (C-iQC-iS) , 

(A7) 

where S is the signal covariance matrix. Its entries are 



(A8) 



where a and (3 are pixels: 1 < a, /3 < A'pix. It could also 
be written as 



(A9) 



where Y^m is a vector of length iVpix containing the val- 
ues of the spherical harmonic Y^m at each pixel. With 
some algebra the expectation value collapses down to: 



(AlO) 



where the window function Wn is: 
1 

Wu^:,{F-% Y^C-iQC-iY,,^. (All) 



2. Computation 

Like the other matrix operations with million-pixel 
maps, direct computation of Wn using Equation (jAlip is 
not feasible; it is 0{N^-^). We have therefore resorted to 
Monte Carlo methods, analogous to those used for trace 
estimation, to simultaneously solve for all of the £ and 
TO terms in Equation (jAlip . Define a random vector z 
of length Npi-^ and with entries consisting of independent 
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random numbers ±1 (i.e. probability 1/2 of being 1 and 
1/2 of being —1). Then (zz-^) ~ 1, so we can write: 

W^i = liF~% ^ yLc-^zz^QC-IY,™^ , 

(A12) 

or 

Wu - ^ ^ (y1„C-1z)(yLc-iC,z)^ . 

(A13) 



To do a Monte Carlo evaluation of the average, we can 
take a random vector z, compute the quantities C~^z and 
C^^CiZ. The latter dominates the computation time, as 
it requires one expensive operation for each power 
spectrum bin, but it is also needed in the Monte Carlo 
evaluation of the Fisher matrix Fij and hence comes with 
no added cost. The inner product Y^^u for any pixel- 
space vector u is the spherical harmonic transform of u, 
for which "fast" 0{Np(^) algorithms exist. We use the 
implementati on o f the spherical harmonic transform of 
Hirata et al. [lTl |. 
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